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ABSTRACT 

The streaming instability (SI) provides a promising mechanism for planetesi- 
mal formation because of its ability to concentrate solids into dense clumps. The 
degree of clumping strongly depends on the height-integrated solid to gas mass 
ratio Z in protoplanetary disks (PPDs). In this letter, we show that the mag- 
nitude of the radial pressure gradient (RPG) which drives the SI (characterized 
by g = tjvk/cs, where tjvk is the reduction of Keplerian velocity due to the RPG 
and Cs is the sound speed) also strongly affects clumping. We present local two- 
dimensional hybrid numerical simulations of aerodynamically coupled particles 
and gas in the midplane of PPDs. Magnetic fields and particle self-gravity are 
ignored. We explore three different RPG values appropriate for typical PPDs: 
q = 0.025,0.05 and 0.1. For each q value, we consider four different particle 
size distributions ranging from sub-millimeter to meter sizes and run simulations 
with solid abundance from Z = 0.01 up to Z = 0.07. We find that a small 
RPG strongly promotes particle clumping in that: 1) At fixed particle size dis- 
tribution, the critical solid abundance Zcrit above which particle clumping occurs 
monotonically increases with q; 2) At fixed Z, strong clumping can occur for 
smaller particles when q is smaller. Therefore, we expect planetesimals to form 
preferentially in regions of PPDs with a small RPG. 
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1. Introduction 



Planetesirnals are super-kilorneter sized bodies that are the building blocks of planets 
( 1Safronovlll969t IChiang &: Youdinll2009l ). yet their formation has long been a mystery. Mil- 
l imeter and centimeter sized particles are routinely observed in protoplanetary di sks (PPDs) 
dWilner et al1l2005[ iRodmann et al.lboOGJ iNatta et al l boO?!: Ihommen et al.lboiol). but parti- 



cles w ith larger size seem difficult to form by coagulation (iBlum fc Wurmll2008l : iGiittler et al. 
20091 ). Solids close to meter size further suffer from rapid radi al drift due to the negative ra- 
dial pressure gradient (RPG) in PPDs (jWeidenschillingl 119771 ) . The gravitational instability 



-2- 



(GI) scenario of planetesimal formation f lGoldreich fc Wardlll973l ) also has difficulties, be- 
cause even without an external source of turbulence, the Kelvin-Helmholtz instability (KHI) 
generated from the dusty midplane layer prevents the onset of GI unless the local height- 
integrated solid to gas mass ratio (Z, her eafter referred to as solid abundance) is about an 
order of magnitude above solar metallicity (IWeidenschillinglll980l : ISekiyalll998l : lYoudin fc Shu 
2002h . 



It was found recently that the dra g interaction between solids and g as leads to a powerful 
insta bility (JGoodman fc Pindorll2000l ). This "streaming instability" (SI. lYoudin &: Goodman 
20051 ) res ults in spontaneous pa rticle clumping from the equilibrium state between solids 
and gas (JNakagawa et al.lll986l ). Numerical simulations demonstrate that the non- linear 



saturation of the SI concentrates particles into dense clumps (jjohansen fc YoudinI 120071 ). 
promoting planetesi mal formation by collective particle self-gravity, bypassing the meter- 
size barrier. Indeed, iJohansen et al.l ( 120071 . |2009| ) showed that planetesimals with sizes of a 
few hundred kilometers form rapidly by the SI from centi-me ter to decimeter sized pebbles 
and rocks, consistent with constraints from the asteroid belt (JMorbidelli et al.l 120091 ). 



This letter com plements our previo us work on the dynamics of particle and gas in the 
midplane of PPDs fJBai fc Stond l2010bl . hereafter BSIO). We consider a wide size distri- 
bution of particles ranging from sub-milli meter up to meter s i zes as an approxi mation for 
the outcome of dust coagulation in PPDs f lBirnstiel et al.ll2010l : IZsom et al.ll2010l ). External 
sources of turbulence suc h as magnetorotational instability are ignored, as appro priate for 
the dead zones in PPDs jGammielll996l : Istone et ahlboooi iBai fc Goodman! I2OO9I ). while SI 
is the main source of turbulence due to particle settling. Particle self- gravity is ignored in 
our simulations, as we focus on the precursor of planetesimal formation: particle clumping. 

The SI is powered by the RPG in the gaseous disk. The RPG reduces the gas orbital 
velocity (in the absence of solids) by a fraction rj of the Keplerian velocity vk = ^r, and is 
characterized by 



q = ttivk/cs = rjr/H, 



g 1 



'l\ 



where c^ is the isothermal sound speed, and Hg = Cg/Q is the scale height of the gaseous disk. 
In BSIO, we assumed q = 0.05 throughout, while in reality q depends on the parameters of, 
and location in, the disk. The strength of the SI turbulence scales with the RPG, which in 
turn affects the particle-gas dynamics in the disk midplane. In particular, we show in this 
letter that particle clumping strongly depends on the RPG, which has important implications 
for planetesimal formation. 



Simulations 



We perform two-diniensional (2D) hybr id simulations of gas and solids using the Athena 
code (IStone et al.ll2008l : lBai fc Stonell2010al ). where gas is treated as a hydrodynamical fluid 
(without magnetic field) on an Eulerian grid, and the solids are treated as superparticles. 
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each representing a swarm of real particles. We model a local patch of the PPDs using the 
shearing sheet approximation. The dynamical equations are written in a reference frame 
corotating at the Keplerian frequency Q at fiducial radius r. We assume axisymmetry, 
and the simulations are performed in the radial-vertical [x-z) plane, where ft is along the z 
direction and x points radially outward. The particles are coupled to the gas via aerodynamic 
drag, characterized by the stopping time tstop, with momentum feedback included. The gas is 
assumed to be isothermal, with midplane density Pg. Vertical gravity Qz = —fl'^z is included 
for both particles and the gas. 

The dynamical equations and simulation setup are identical to those in BSIO. Specif- 
ically, we consider a particle size distribution that is discretized into a number of particle 
size bins, with each bin covering half a dex in r^ = fitstop? bounded by the minimum and 
maximum stopping time Tmin and Tmax- In what follows, we label our simulations with names 
of the form Kmn, where m, n are integers obtained by Tmin = 10~™ and Tmax = lO"", thus 
run Kmn uses 2{m — n) + 1 particle size bins (or particle species). We use a variety of grid 
resolution and box sizes to capture the fastest growing modes of the SI, and we use 10^ 
particles per species in all our simulations. We assume uniform mass distribution across all 
particle size bins, with total solid abundance Z. As in BSIO, we consider four groups of 
runs, R41, R21, R30 and RIO. 

For each particle size distribution, we perform simulations with three different values of 
the RPG parameter: q = 0.025,0.05 and 0.10. According to equations (4) and (7) in BSIO, 
the dependence of q on disk parameters such as temperature and mass is relatively weak, 
thus our range of q covers a large parameter space of disk models. For each set of runs, we 
perform a series of simulations with different Z values, starting at Z = 0.01 and increasing 
the value by 0.01 for each new run, until strong particle clumping occurs oi Z = 0.07. Our 
simulation run parameters are summarized in Table [U They are identical to the 2D run 
parameters in BSIO, but use different values for q and a larger range in Z. Since the natural 
length scale of the SI is rjr, we use smaller (bigger) simulation box sizes for smaller (larger) 
q values so that rjr is resolved by an equal number of grid cells in each series of runs. All 
simulations are run for at least 900^2"^. 

Three-dimensional (3D) s imulation with the inclusion of the azimuthal dimension is 



necessary to capture the KHI ( 1Chiangil2008t iBarrancd 120091 ). which is mainly caused by the 



vertical shear in the gas azimuthal velocity. Our simulations are 2D rather than 3D for two 
reasons. First, for our adopted particle size distribution, the turbulence generated from the 
SI stops particles from setthng before the onset KHI (BSIO), at least for the q = 0.05 case. 



Second , in order to properly resolve the SI, relatively high resolution is required ( jBai fc Stone 



2010al ). and particle clumping does depend on resolution|^. Simulations in 3D with such high 



resolution are too costly due to the small box size (thus time step), especially for q = 0.025. 



^For example, Run R30Z3-3D in BSIO has no particle clumping in our standard resolution, but shows 
clumping at lower resolution. 
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3. Particle Clumping 



The simulations saturate in about 50-100 orbits. Particle settling triggers the SI, and 
particles with different stopping times are maintained at different heights determined by the 
balance between settling and turbulent diffusion. Most interestingly, SI effi ciently concen- 



trate particles into dense clumps when the solid abundance is sufficiently high ( iJohansen et al 



20091 . BSIO). In Figured] we plot the evolution of maximum particle density Pp,max for the 
R30 runs with different values of q and Z. There is a clear dichotomy in the evolution: either 
the maximum particle density stays at a relatively small value (pp,max ^ 50pg), or it reaches 
as high as lO^p^, indicative of strong clumping. A particle clump becomes gravitationally 
bound when its density exceeds the Roche density (see equation (18) of BSIO), which also is 
of the order lO^pg for typical PPDs. Therefore, particle clumping is a prelude to planetesimal 
formation via gravitational collapse. For each value of g, the transition from non-clumping to 
clumping is sharp as Z gradually increases, and we can define a critical solid abunance Zcrit, 
where strong particle clumping occurs for Z > ^crit- Comparing different panels indicates 
that Zcrit monotonically increases with q. 

To better demonstrate the effect of the disk RPG on planetesimal formation, we plot the 
maximum particle density as a function of solid abundance Z for each group of simulations 
with all values of q in Figure [21 The maximum density is taken from the largest value of 
Pp,max over the last 20 orbits of the simulations, when all the runs are fully saturated. We 
take pp = 500pg (which is the same order as the Roche density) as a rough indicator of 
planetesimal formation, as shown by the dotted lines. For each value of q, one can determine 
Zcrit at the intersection of the dashed and corresponding solid lines. We see that for all 
four groups of runs, Z^rit monotonically increases with q. Moreover, Zcnt depends on q more 
sensitively when the particle size is on average smaller. When q = 0.025, Z^nt is about 0.015 
for all our four groups of runs. For g = 0.1, run R41 does not show particle clumping below 
Z = 0.07, for runs R21 and R30 Zcnt is about 0.06, while for runs RIO Zcnt drops to below 
0.05. 

Particle clumping is a highly non-linear effect due to the SI, and is more likely to develop 
when the average disk midplane solid to gas mass ratio e is large. If one assumes D to be the 
midplane vertical diffusion coefficient of the SI turbulence, particles with stopping time Tg 
would settle to a layer with thickness of the order Hp ^ y^D/Q,Ts. The diffusion coefficient 
D depends on the particle size distribution, e ~ ZHg/Hp and q. Without vertical gravity, 
the only length scale in the problem is rjr = qHg, thus D oc q^. If one assumes D oc e", we 
obtain by generalizing the toy model in BSIO (see their equation (16) and (17)) 

Hp OC Z"/(°+2)^2/(a+2) _ ^2) 

Since one expects a < for relatively large e > 1, we see that Hp not only depends sensitively 
on Z, as shown in BSIO, it depends even more sensitively on q. This critical dependence 
makes the average particle to gas mass ratio e at midplane quickly increases as q decreases, 
promoting strong particle clumping, which explains the trend in Figure [2l The more sensitive 
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Table 1: Run parameters. 



Run 



R41 



R21 



R30 



RIO 



Q. 



lOOZ r, 



min 



r, 



max 



L/^ X L/r. 



N.y<N, 



0.025 1..3 
0.05 1..5 
0.1 1..7 



10" 



10 



0.05 X 0.15 256 x 768 
-1 0.1x0.3 
0.2 X 0.4 



256 X 768 
256 X 512 



0.025 1..3 
0.05 1..4 
0.1 1..7 



10 



-2 



10' 



0.05 X 0.15 
0.1 X 0.3 
0.2 X 0.4 



0.025 1..3 
0.05 1..3 
0.1 1..7 



10' 



0.2 X 0.3 
0.2 X 0.3 
0.2 X 0.3 



0.025 1..3 
0.05 1..3 
0.1 1..6 



10' 



0.2 X 0.3 
0.2 X 0.3 
0.2 X 0.3 



^ Domain size, in unit of gas scale height Hg = Cg/Q. 
^ Grid resolution. 



256 X 768 
256 X 768 
256 X 512 



256 X 512 
256 X 384 
256 X 384 



256 X 512 
256 X 384 
256 X 384 




T . =10" 

mm 



T =1 
max 



-Z=0.01 
-Z=0.02 
-Z=0.03 
-Z=0.04 
Z=0.05 
Z=0.06 
-Z=0.07 



200 400 600 800 200 400 600 800 200 



t (ii"^) 



t(n 



-1\ 
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t(n" 



600 



800 



Fig. 1. — Time history of the maximum particle density in simulations with a range of 
particle sizes from r^ = 10"^ to one (runs R30). Left, middle and right panels show the 
results from q = 0.025, 0.05 and 0.1 respectively. Runs with different solid abundance Z are 
labeled with different colors. 
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dependence of Zcrit on q for smaller particles can be interpreted as a result of the power-law 
index a tending to be more negative for smaller r^. 

Our 2D simulations do not contain the azimuthal dimension, which suppresses KHI. 
Nevertheless, we can check whether KHI would occur before particles settle into a sufficiently 
thin layer to trigger strong SI by plotting the vertical profile of the Richardson number for 
simulations at critical solid abundance for all q values, see Figure El The Richardson number 
is calculated by the method described in §3.1 of BSIO. Alt hough the Richardson nu mber 



alone does n ot deterniine t he stability when Coriolis force (JGomez fc Ostrikerll2005l ). and 



radial shear (JBarrancdl2009l ) are pre sent, one may still take Rirrn. = 0.1 as an approximation 



for the critical value for instability (jChiang|2008l : iLee et al.l 120101 ). We see in Figure [3] that 



for most of these runs, the Richardson number across the disk is above the critical value, 
meaning that the system is KH stable at ^crit- Moreover, as q increases, the Richardson 
number at Zcrh decreases. In particular, for R30 and RIO runs with g = 0.1, the Richardson 
number at the disk midplane drops below 0.1, which may be subject to KHI. Therefore, the 
condition for strong clumping for these runs may be more stringent than shown in Figure [2l 



4. Discussions and Conclusions 

In this letter, we have demonstrated that a small RPG strongly favors particle clumping 
and planetesimal formation via three effects. First, for a fixed particle size distribution, the 
critical solid abundance Zcnt above which strong clumping occurs is reduced. Second, at 
fixed solid abundance, strong clumping occurs for smaller particles. Third, KHI is less likely 
to be triggered (and therefore will not suppress particle clumping) at Z = Zcrit- Moreover, 
a smaller RPG also implies smaller radial drift and collision velocities, promoting grain 
growth, strengthening our second point above. A small RPG also favors the GI scenario of 
planetesimal formation (relevant when all particles are strongly coupled to the gas) because 
there is less free energy available in the vertical shear □. 

One caveat from our 2D simulations is that the condition for particle clumping is slightly 
more stringent in 3D than in 2D (BSIO). Therefore, the values of Zcnt obtained in this letter 
may be considered as a lower bound. Nevertheless, the conclusions in this letter are robust. 
In fact, Johansen et al. (2007, see their supplemental material) also found that planetesimal 
formation is eased with smaller pressure gradient using 3D simulations. 

Combined with the results in BSIO, we have shown that strong particle clumping is 
favored when there is: 1) a small RPG; 2) large solid abundance; 3) large solid size; and 4) a 
dead zone. These results indicate that planetesimals preferentially form in specific locations 
in PPDs. The global structure and evolution of PPDs is of crucial importance: it deter- 
mines the RPG profile, the radial transport of particles which leads to enhancement of solid 



^We expect the threshold abundance for GI to operate is larger than Z^rit- 
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Fig. 2. — The maximum particle density that can be achieved by SI from all our simulation 
runs (for R41, R21, R30, RIO from left to right panels). In each panel, we plot the maximum 
particle density as a function of solid abundance Z for each value of q (black solid, blue 
dashed and red dash-dotted for q = 0.025,0.5 and 0.1 respectively). The dotted line at 
Pp,max = 500pg indicates the adopted threshold for strong particle clumping, and we use 
larger symbols to indicate the first runs that show strong clumping as Z increases. 
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Fig. 3. — The Richardson number profile for all simulation runs (for R41, R21, R30, RIO 
from left to right panels) with solid abundance just above Zcnt (the runs with enlarged 
symbols in Figure |2]). Results from q = 0.025,0.05 and 0.1 are shown in black solid, blue 
dashed and red dash-dotted curves respectively. The dotted line indicated an approximate 
estimate of the critical Richardson number Ricrit = 0.1. 



abundance in the inner disk fJYoudin fc Shull2002l ). and grain coagulation that determines 
the particle size distribution. R ecent models on the structure and evolution of PPDs (e.g., 
Jin fc Suill2010l : IZhu et al.ll2010l ) generally show much more complicated structures than the 
simple MMSN or a-disk models due to non-steady state accretion as well as the presence 
of dead zones. In ad dition, RPG may reach zero in local pressure bumps at the sii ow line 
(JKretke fc Liru 120071 ) or the inner edge of the dead zone (JDzyurkevich et al.l |2010| ) which 
can also be the preferred site for planetesimal formation. These results all-together provide 
useful constrairits on the in itial conditions for the formation of planet embryos and planets 
kokubo fc Idall2000l . [2002l ). 
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